///\todo upgrade

//////////////////////////////////////////////////////////////////////////////////
//                                                                              //
//  This file is part of the buola project (https://code.google.com/p/buola/).  //
//                                                                              //
//  Copyright(c) 2007-2012 Xavi Gratal                                          //
//  gratal AT gmail DOT com                                                     //
//                                                                              //
//  Buola is free software: you can redistribute it and/or modify               //
//  it under the terms of the GNU General Public License as published by        //
//  the Free Software Foundation, either version 3 of the License, or           //
//  (at your option) any later version.                                         //
//                                                                              //
//  Buola is distributed in the hope that it will be useful,                    //
//  but WITHOUT ANY WARRANTY; without even the implied warranty of              //
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               //
//  GNU General Public License for more details.                                //
//                                                                              //
//  You should have received a copy of the GNU General Public License           //
//  along with buola.  If not, see <http://www.gnu.org/licenses/>.              //
//                                                                              //
//////////////////////////////////////////////////////////////////////////////////

#ifndef _BUOLA_MAT_CHOLESKY_H_
#define _BUOLA_MAT_CHOLESKY_H_

#include <buola/mat.h>
#include <buola/mat/detail/eigen_wrapper.h>
#include <buola/mat/views.h>
#include <Eigen/Cholesky>

namespace buola { namespace mat {

DECLARE_NAME(FromL)
    
template<typename tScalar>
class CLLT
{
public:
    CLLT()=default;
    CLLT(const CLLT&)=default;
    CLLT(CLLT&&)=default;
    CLLT& operator=(const CLLT&)=default;
    CLLT& operator=(CLLT&&)=default;
    
    template<typename tExpr>
    CLLT(NFromL,tExpr &&pExpr)
        :   mL(std::forward<tExpr>(pExpr))
    {}
    
    detail::CTriangularView<const mat::CMat<tScalar>&,NLower> L() const    {   return {mL};  }
    detail::CTriangularView<detail::CTranspose<const mat::CMat<tScalar>&>,NUpper> U() const  {   return {nConstructSub,mL};  }
    
    template<typename tDerived>
    detail::OMatrixType<tDerived> Solve(const CExpr<tDerived> &pExpr)
    {
        return U().Solve(L().Solve(pExpr.Derived()));
    }
    
private:
    CMat<tScalar> mL;
};
    
template<typename tDerived>
CLLT<detail::OScalar<tDerived>> cholesky(const CExpr<tDerived> &pExpr)
{
    const auto &lE=pExpr.Derived();
    assert(lE.Rows()==lE.Cols());
    mat::CMat<detail::OScalar<tDerived>> lL(lE.Rows(),lE.Rows());
    detail::eigen_map(lL)=detail::eigen_map(detail::eigen_mappable(lE)).llt().matrixL();
    return {nFromL,lL};
}

/*namespace mat*/ } /*namespace buola*/ }

#endif
